"""
多变量lstm预测，还没有修改，不可用
"""
from math import sqrt
from numpy import concatenate
from matplotlib import pyplot
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from pandas import Series
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import mean_squared_error
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from pandas import datetime
import numpy as np
import matplotlib.pyplot as plt
from keras.callbacks import EarlyStopping
from keras.layers import Bidirectional
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score


def parse(x):
    return datetime.strptime(x, '%Y %m')


# convert series to supervised learning
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
    n_vars = 1 if type(data) is list else data.shape[1]
    df = DataFrame(data)
    cols, names = list(), list()
    # input sequence (t-n, ... t-1)
    for i in range(n_in, 0, -1):
        cols.append(df.shift(i))
        names += [('var%d(t-%d)' % (j + 1, i)) for j in range(n_vars)]
    # forecast sequence (t, t+1, ... t+n)
    for i in range(0, n_out):
        cols.append(df.shift(-i))
        if i == 0:
            names += [('var%d(t)' % (j + 1)) for j in range(n_vars)]
        else:
            names += [('var%d(t+%d)' % (j + 1, i)) for j in range(n_vars)]
    # put it all together
    agg = concat(cols, axis=1)
    agg.columns = names
    # drop rows with NaN values
    if dropnan:
        agg.dropna(inplace=True)
    return agg


def mean_absolute_percentage_error(y_true, y_pred):
    # 平均绝对百分比误差（MAPE）的计算
    y_true, y_pred = np.array(y_true), np.array(y_pred)
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100


# evaluate the RMSE for each forecast time step
def evaluate_forecasts(test, forecasts):

    rmse = sqrt(mean_squared_error(test, forecasts))
    mae = mean_absolute_error(test, forecasts)
    mape = mean_absolute_percentage_error(test, forecasts)
    r2 = r2_score(test, forecasts)
    print(' RMSE: %f' % rmse)
    print(' MAE: %f' % mae)
    print(' MAPE: %f' % mape)
    print(' R2: %f' % r2)

    return r2


if __name__ == '__main__':

    # 加载数据
    dataset = read_csv(
        'E:\lyf_ML_Drought\coding\ML_Drought_Prediction\indices_caculate\\result\multi_spei_csv\SPEI-12\Multi_SPEI-12_52533.txt',
        header=0, parse_dates=[['year', 'month']], index_col=0, date_parser=parse)
    dataset.index.name = 'time'
    dataset.drop(columns=['average_air_pressure', 'average_water_air_pressure'], inplace=True)
    print(dataset)
    values = dataset.values

    # 将字符数据编码为数字数据，这里用不到
    # encoder = LabelEncoder()
    # values[:,4] = encoder.fit_transform(values[:,4])

    # 确保数据都是浮点型
    values = values.astype('float32')

    # 数据缩放
    scaler = MinMaxScaler(feature_range=(-1, 1))
    scaled = scaler.fit_transform(values)

    # 将数据转化为监督数据,再把不需要的列删除
    reframed = series_to_supervised(scaled, 1, 1)
    reframed.drop(reframed.columns[[8, 9, 10, 11, 12, 13]], axis=1, inplace=True)
    # print(reframed)

    # split into train and test sets
    values = reframed.values
    train_step = int(0.9*len(values))
    train = values[:train_step, :]
    test = values[train_step:, :]

    # split into input and outputs
    train_X, train_y = train[:, :-1], train[:, -1]
    test_X, test_y = test[:, :-1], test[:, -1]
    # print("train_X",train_X)
    # print("train_y",train_y)
    # print("test_X",test_X)
    # print("test_y", test_y)

    # reshape input to be 3D [samples, timesteps, features]
    train_X = train_X.reshape((train_X.shape[0], 1, train_X.shape[1]))
    test_X = test_X.reshape((test_X.shape[0], 1, test_X.shape[1]))
    # print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)

    # design network
    model = Sequential()
    model.add(LSTM(50, input_shape=(train_X.shape[1], train_X.shape[2]), stateful=False))
    model.add(Bidirectional(LSTM(7, activation='tanh'), input_shape=(train_X.shape[1], train_X.shape[2])))
    model.add(Dense(1))
    model.compile(loss='mae', optimizer='adam')
    # fit network
    # early_stop = EarlyStopping(monitor='loss', patience=10, verbose=1)
    # history = model.fit(train_X, train_y, epochs=3000, batch_size=1, validation_data=(test_X, test_y), verbose=2, shuffle=False,, callbacks=[early_stop])
    history = model.fit(train_X, train_y, epochs=3000, batch_size=50, verbose=2,
                        shuffle=False)
    # plot history
    # pyplot.plot(history.history['loss'], label='train')
    # pyplot.plot(history.history['val_loss'], label='test')
    # pyplot.legend()
    # pyplot.show()

    # make a prediction
    yhat = model.predict(test_X)
    test_X = test_X.reshape((test_X.shape[0], test_X.shape[2]))

    # invert scaling for forecast
    inv_yhat = concatenate((yhat, test_X[:, 1:]), axis=1)
    inv_yhat = scaler.inverse_transform(inv_yhat)
    inv_yhat = inv_yhat[:, 0]

    # invert scaling for actual
    test_y = test_y.reshape((len(test_y), 1))
    inv_y = concatenate((test_y, test_X[:, 1:]), axis=1)
    inv_y = scaler.inverse_transform(inv_y)
    inv_y = inv_y[:, 0]

    # evaluate forecasts
    r2 = evaluate_forecasts(test_y, yhat)

    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(range(len(test_y)), test_y, c="black")
    ax.plot(range(len(yhat)), yhat, c="red")
    # x_ticks = ax.set_xticks(range(len(true_list) / 24, len(true_list), len(true_list) / 12))
    # x_labels = ax.set_xticklabels(range(1, 13, 1), fontsize="small")
    ax.set_xlabel("Month")
    ax.set_ylabel("SPEI")
    ax.legend(['true', 'pred'])
    plt.show()
